
#ifndef AFTERMATH_PROBABILITY_STATISTICS_BUILDER_HPP_INCLUDED
#define AFTERMATH_PROBABILITY_STATISTICS_BUILDER_HPP_INCLUDED

#include "../type_concepts.hpp"

#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <type_traits>
#include <utility>

namespace aftermath
{
    namespace probability
    {
        namespace detail
        {
            template <typename t_key_type, typename t_probability_type,
                bool t_is_enabled = type_impl::has_less<t_key_type>::value>
            struct order_statistic
            {
                static const bool is_enabled = t_is_enabled;

            private:
                bool m_initialized = false;
                t_key_type m_min, m_max;

            public:
                void clear()
                {
                    this->m_initialized = false;
                }

                void observe(const t_key_type& key, bool& is_new_min, bool& is_new_max)
                {
                    if (!this->m_initialized)
                    {
                        this->m_min = key;
                        this->m_max = key;
                        is_new_min = true;
                        is_new_max = true;
                        this->m_initialized = true;
                    }
                    else if (key < this->m_min)
                    {
                        this->m_min = key;
                        is_new_min = true;
                    }
                    else if (this->m_max < key)
                    {
                        this->m_max = key;
                        is_new_max = true;
                    }
                }

                const t_key_type& get_min() const
                {
                    return this->m_min;
                }

                const t_key_type& get_max() const
                {
                    return this->m_max;
                }

                template <typename t_dictionary_type>
                t_probability_type cdf(const t_key_type& key, typename t_dictionary_type::mapped_type norm, const t_dictionary_type& data) const
                {
                    if (key < this->m_min) return 0.0;
                    if (!(key < this->m_max)) return 1.0;

                    typedef typename t_dictionary_type::mapped_type count_type;

                    count_type cumulative_count = count_type(); // zero
                    for (const auto& item : data)
                    {
                        if (key < item.first) break;
                        cumulative_count += item.second;
                    }
                    return cumulative_count / static_cast<t_probability_type>(norm);
                }

                template <typename t_dictionary_type>
                const t_key_type& percentile(t_probability_type probability, typename t_dictionary_type::mapped_type norm, const t_dictionary_type& data) const
                {
                    //if (probability < 0.0) cresult = std::error_condition(-1, std::out_of_range);
                    //if (probability > 1.0) cresult = std::error_condition(1, std::out_of_range);  
                    //cresult.clear();

                    if (probability <= 0.0) return this->m_min;
                    if (probability >= 1.0) return this->m_max;

                    typedef typename t_dictionary_type::mapped_type count_type;

                    count_type p = count_type(); // zero
                    // Think in counts rather than probabilities.
                    probability *= norm;
                    count_type threshold = static_cast<count_type>(probability);
                    if (std::is_integral<count_type>::value && threshold < probability) threshold++; // For integral types: we need a ceiling, not floor.

                    for (const auto& item : data)
                    {
                        p += item.second;
                        if (p >= threshold) return item.first;
                    }
                    return this->m_max;
                }
            };

            template <typename t_key_type, typename t_probability_type>
            struct order_statistic<t_key_type, t_probability_type, false>
            {
                static const bool is_enabled = false;

                void clear() const { }
                void observe(const t_key_type&, bool&, bool&) const { }
            };

            template <typename t_key_type, typename t_count_type, typename t_sum_type, typename t_mean_type,
                bool t_is_enabled = 
                    type_impl::has_product<t_key_type, t_count_type>::value &&
                    type_impl::has_sum<t_sum_type, t_key_type>::value &&
                    type_impl::has_divide<t_mean_type, t_count_type>::value>
            struct linear_statistic
            {
                static const bool is_enabled = t_is_enabled;

            private:
                t_sum_type m_sum = t_sum_type();

            public:
                void clear()
                {
                    this->m_sum = t_sum_type();
                }

                void observe(const t_key_type& key, t_count_type repeat)
                {
                    this->m_sum += static_cast<t_sum_type>(key * repeat);
                }

                const t_sum_type& get_sum() const
                {
                    return this->m_sum;
                }

                t_mean_type get_mean(t_count_type norm) const
                {
                    return static_cast<t_mean_type>(this->m_sum) / norm;
                }

                //// NOTE: here comes the quadratic part of a linear statistic *facepalm*
                //template <typename t_dictionary_type>
                //t_mean_type compute_variance(t_count_type norm, const t_dictionary_type& data) const
                //{
                //    return variance_statistic<t_key_type, t_count_type, t_sum_type, t_mean_type, is_quadratic>::compute_variance(norm, data, this->get_mean(n));
                //}

                //template <typename t_dictionary_type>
                //t_mean_type compute_standard_deviation(t_count_type norm, const t_dictionary_type& data) const
                //{
                //    return variance_statistic<t_key_type, t_count_type, t_sum_type, t_mean_type, is_quadratic>::compute_standard_deviation(norm, data, this->get_mean(n));
                //}
            };

            template <typename t_key_type, typename t_count_type, typename t_sum_type, typename t_mean_type>
            struct linear_statistic<t_key_type, t_count_type, t_sum_type, t_mean_type, false>
            {
                static const bool is_enabled = false;
                
                void clear() const { }
                void observe(const t_key_type&, std::size_t) const { }
            };

            template <typename t_key_type, typename t_count_type, typename t_sum_type, typename t_mean_type,
                bool t_is_enabled =
                    linear_statistic<t_key_type, t_count_type, t_sum_type, t_mean_type>::is_enabled &&
                    type_impl::has_difference<t_key_type, t_mean_type>::value &&
                    std::is_scalar<typename type_impl::get_product_result<t_key_type, t_key_type>::type>::value>
            struct variance_statistic
            {
                static const bool is_enabled = t_is_enabled;

                template <typename t_dictionary_type>
                static double compute_variance(t_count_type norm, const t_dictionary_type& data, const t_mean_type& mean)
                {
                    if (norm == 0) return std::numeric_limits<double>::signaling_NaN();

                    double variance_sum = t_mean_type(); // zero.
                    for (const auto& item : data)
                    {
                        variance_sum += ((item.first - mean) * (item.first - mean)) * item.second;
                    }
                    return variance_sum / norm;
                }

                template <typename t_dictionary_type>
                static double compute_standard_deviation(t_count_type norm, const t_dictionary_type& data, const t_mean_type& mean)
                {
                    return std::sqrt(compute_variance(norm, data, mean));
                }
            };

            template <typename t_key_type, typename t_count_type, typename t_sum_type, typename t_mean_type>
            struct variance_statistic<t_key_type, t_count_type, t_sum_type, t_mean_type, false>
            {
                static const bool is_enabled = false;
            };
        }
    }
}

#endif // AFTERMATH_PROBABILITY_STATISTICS_BUILDER_HPP_INCLUDED
